OPEN 3 ACCESS Freely available online 



•0-PLOS I o-^E 



Modeling the Potential Impact of Host Population (Jfi 
Survival on the Evolution of M, tuberculosis Latency cros^Mark 

Nibiao Zheng\ Christopher C. Whalen^, Andreas HandeP* 

1 Institute of Bioinformatics, University of Georgia, Athens, Georgia, United States of America, 2 Department of Epidemiology and Biostatistics, College of Public Health, 
University of Georgia, Athens, Georgia, United States of America 

Abstract 

Tuberculosis (TB) is an infectious disease with a peculiar feature: Upon infection with the causative agent, Mycobacterium 
Tuberculosis (MTB), most hosts enter a latent state during which no transmission of MTB to new hosts occurs. Only a fraction 
of latently infected hosts develop TB disease and can potentially infect new hosts. At first glance, this seems like a waste of 
transmission potential and therefore an evolutionary suboptimal strategy for MTB. It might be that the human immune 
response keeps MTB in check in most hosts, thereby preventing it from achieving its evolutionary optimum. Another 
possible explanation is that long latency and progression to disease in only a fraction of hosts are evolutionary beneficial to 
MTB by allowing it to persist better in small host populations. Given that MTB has co-evolved with human hosts for millenia 
or longer, it likely encountered small host populations for a large share of its evolutionary history and had to evolve 
strategies of persistence. Here, we use a mathematical model to show that indeed, MTB persistence is optimal for an 
intermediate duration of latency and level of activation. The predicted optimal level of activation is above the observed 
value, suggesting that human co-evolution has lead to host immunity, which keeps MTB below its evolutionary optimum. 
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Introduction 

Tuberculosis (TB) is an infectious disease caused by the 
bacterium Mycobacterium tuberculosis (MTB). MTB lias been 
infecting humans for a long time, at least millenia and likely even 
longer [1]. As such, MTB had to evolve evolutionary strategies 
that allowed it to persist in small groups of human hosts. An 
interesting question is if one can see signatures of this potential 
evolutionary adaptation to small groups of human hosts in MTB's 
"life history" today. 

Upon infection with MTB, most hosts enter a latent state. Those 
hosts do not show signs of disease but do harbor MTB, which can 
activate and lead to disease at any future time [2-4] . It might take 
a long time before activation occurs, and the majority of TB 
infected hosts die from other causes besides TB without ever 
developing disease [5]. Hosts latently infected with MTB cannot 
infect others. Therefore, at first glance, latency does not seem to be 
beneficial for MTB. One possible explanation for the long latency 
and the fact that activation to disease only occurs in a fraction of 
hosts is based on the human immune response [6-8]. It is known 
that a competent immune response is needed to contain infection 
and avoid disease, as dramatically illustrated in HIV infected hosts 
with weakened immune responses, who activate at much higher 
rates [9-11]. This suggests that the co-evolutionary arms race 
between MTB and its human hosts has lead to a host immune 
response that can successfully contain MTB in a state of 
suboptimal fitness [12,13]. The fact that there seems to be local 
adaptation on both the pathogen and host side lends support to 
this idea [13-17]. 



Although host immune pressure is a plausible explanation for 
reduced activation rates, there are also arguments against this idea. 
Pathogens like MTB replicate rapidly (compared to their human 
hosts) and often reach large population sizes; both features are 
known to foster rapid evolution, especially under strong selective 
pressure [18,19]. This is at vivid display in the evolution of drug 
resistance for MTB and many other pathogens [20,2 1] . Therefore, 
one could argue that if long latency and low activation rates were 
evolutionary strongly suboptimal for MTB, evolution would have 
led to higher rates of activation. Instead, as has been suggested 
previously [22,23], long latency and low activation rates might be 
strategies that are evolutionary beneficial to MTB by increasing its 
fitness. 

There is no single way to quantify the fitness of an organism. 
For infectious diseases, an often used measure of fitness is 
transmissibHity, usually defined by the reproductive number, 
[24,25]. It can be shown that in direct competition of two strains 
(under equilibrium conditions), the strain with the higher Rq 
outcompetes the one with the lower [25]. However, in the 
absence of direct competition, strains with lower R^ might 
sometimes be more advantaged as they can better avoid local 
extinction and therefore win through indirect competition against 
strains with higher transmissibUity [26-29]. Therefore, an 
alternative way to capture fitness of pathogens is by quantifying 
their ability to persist in a host population. 

The impact of transmissibUity and persistence on overall fitness 
has been an active area of research [30-34] . The importance of 
different measures of fitness such as these depends on the situation. 
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Figure 1. Flow diagram for the TB transmission model. Susceptible hosts (S) are born at rate X with birth rate saturating at high population 
density, W^- Hosts in each compartment die at a background mortality rate m„. Infection of susceptibles occurs at rate b upon contact with TB 
infectious hosts (I). After TB infection, a fraction f of hosts develop active TB in a short time (fast progression), the remainder enter the latent stage. 
Latent hosts (L) can activate and develop TB disease sometime later at rate a (slow progression). Latent hosts can also develop active TB after re- 
infection with TB, with probability of this occurring compared to fast progression of susceptibles reduced by a factor k. Infectious, diseased hosts (I) 
might regress to the latent stage at rate w. Hosts with active disease die due to diseased induced mortality at rate m^. 
doi:1 0.1 371/journal.pone.0105721 .gOOl 



Since MTB is an ancient organism that iias infected humans for 
mUlenia [12,35-38], it had to evolve strategies that allowed it to 
persist in relatively small, spatially structured host populations. 
Under such circumstances, fitness benefits due to non-extinction 
might have had an important impact on overall evolutionary 
dynamics [28,39,40]. 

It has been previously proposed that a prolonged latency might 
have been one persistence strategy [23,41]. Here, we use a 
mathematical model to explore this idea. We introduce a measure 
of persistence of MTB in a population and investigate how well 
MTB can persist as a function of the latent period. We find that 
TB persistence is optimal for an intermediate duration of latency 
and level of activation. We also find that the optimal level of 
activation is above the observed value, suggesting that host 
immunity plays some role in keeping MTB below its optimal level 
of fitness. 

Methods 

Mathematical model description 

We use a compartmental mathematical model formulated as a 
set of ordinary differential equations to describe the population- 
level infection dynamics of TB. The model is shown schematically 
in figure 1. Our model is similar to other recently studied TB 
models [42,43]. We consider 3 types of hosts (compartments) in 
our model: susceptible hosts, S, latently infected hosts, L, that 
harbor MTB but are not infectious and show no signs of disease, 
and infectious hosts with active disease, /. We keep the model 



simple and do not distinguish between features such as age-related 
differences (e.g. children versus adults). Such additional details 
could be included in fiarther more detailed models. The total 
population size is N — S+L+I. New hosts enter the system at a 
maximum rate 1 per person, this rate saturates once the 
population reaches some maximum level, N„^. All hosts die due 
to causes other than TB after an average lifespan of years. 
Uninfected hosts can become infected through contact with an 
infectious host at rate b. The infection process is modeled in a 
density dependent manner [43,44]. After infection, a small 
fraction / of hosts rapidly develop disease (fast progression) 
[1 1,45,46] and move into the active disease compartment, /. The 
majority of hosts enter the latent state, L (slow progression). 
Latently infected hosts can convert to infectious, diseased hosts 
later in their life at rate a (activation) or through reinfection, with 
the chance of disease due to reinfection reduced by a factor k [47- 
49]. Infectious, diseased hosts either die due to disease after on 
average l/ma years, or regress and return to the latent stage at rate 
w. Following previous models, we assume that recovered 
individuals do not fuUy clear the infection but instead return to 
the latent stage [43,50-52]. Table 1 summarizes the model 
variables and parameters and provides references for the 
parameter values used for most of our analysis. The model 
equations are 



dS 
dt 



--XN(\- 



N 
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— = (1 -f)bIS + wl -aL- kfhIL - rrinL (2) 



— =fhIS + aL + kfbIL-wI-(m„+md)I (3) 

We implemented all model simulations in R [53]. The code to 
reproduce the simulations described here is available from the 
author's webpage at http:/ /handelgroup. uga.edu/resources. htm. 

Results 

Persistence of MTB as a function of latency 

Our main outcome of interest is the potential of MTB to persist 
in a host population and how persistence potential depends on the 
duration of the latent period and fraction of latent TB hosts that 
activate. While pathogens that have an environmental stage as an 
important component of their transmission cycle can persist even 
in the absence of any infected hosts [54-58], for pathogens where 
direct transmission is the main important component, persistence 
(non-extinction) requires the continued presence of infected hosts. 
For MTB, extinction occurs if no more latendy infected hosts, L, 
and infectious, diseased hosts, /, are present. Persistence is more 
likely (i.e. extinction is less likely) as the number of infected hosts in 
the population increases [59-61]. Since for MTB, only a fraction 
of latently infected hosts wiU become infectious and contribute to 
transmission, a better measure of persistence is given by the 
quantity P, defined as P = I + a.L, where / and L are the number 
of infectious and latently infected hosts, and a represents the 
fraction of latentiy infected hosts that wiU develop TB disease, 
become infectious and are able to transmit. For our model, a is 
given by the ratio of the rate of latently infected hosts that move on 
to develop active TB, a + kjbl, divided by the total rate at which 
latendy infected hosts leave the latent stage, a+kfbl+m„, i.e. 
a = (a + kfil) /(fl + kjbl + in„). 



Table 1. Initial conditions and parameter values. 



Persistence as defined by P, especially at the steady state, 
provides a useful measure for the ability of MTB to not go extinct 
in the population (see our comparison with a stochastic model 
below). We run simulations of our model for different values of the 
activation rate, a, and record / and L at steady state and use this to 
compute P as function of a. Figure 2 A shows that optimal 
persistence is achieved at intermediate rates of activation. 
Figure 2B shows individually the three components that make 
up P. The number of infectious hosts at steady state, as well as the 
fraction of hosts activating, cc, increase with increasing activation 
rate. The number of latent hosts first increases and then decreases 
with larger activation rate. The combination of these three 
quantities leads to a maximum for persistence P at intermediate 
values. 

Optimal persistence at an intermediate rate of activation also 
implies that instead of having every latent host activate and 
become infectious, it is beneficial for the pathogen to let some 
infections "go to waste" by way of latent hosts dying before they 
become infectious. This helps with overall persistence and is worth 
the "loss" of a fraction of latent hosts due to natural death before 
they activate and are able to transmit. Figure 2C shows persistence 
as a function of the fraction of hosts that activate, a. The figure 
also illustrates another interesting point: The optimal fraction of 
hosts that activate is slightly above 50% given the chosen model 
parameters. This is higher than the —10-20% observed [62,63], 
suggesting that MTB is not able to achieve the activation rate that 
would optimize its persistence. This might be attributable to the 
host immune response playing a role at reducing activation. 

The impact of parameter value uncertainty on optimal 
persistence and activation 

To investigate how sensitive results reported in the previous 
section are to changes in parameter values, we performed an 
uncertainty analysis and sampled the model parameters using 
Latin Hypercube Sampling [64-67]. For each parameter, we 
considered uniform distributions ranging between 0.5 and 2 times 
the base parameter value shown in table 1. For each parameter 
sample we computed the values Po, flo and o(o, i.e. the optimal level 
of persistence and the activation rate and fraction at which the 





Symbol 


Interpretation 


Value 


Source and Comments 




maximum population size 


1 


population size normalized to 1 


So 


initial susceptible hosts 


0.6 


calculated as (\—m,JX)N,^,, to obtain steady population size in the 
absence of TB 


Lo 


initial latent hosts 


0 


arbitrary choice 


lo 


initial infectious, diseased hosts 


So/1000 


one infected per 1,000 hosts 


X 


maximum birth rate 


0.05 per year 


50 births per 1,000 hosts, representing a high birth rate scenario 


m„ 


natural mortality rate 


0.02 per year 


assuming a life-span of 50 years for healthy hosts 




disease-induced mortality rate 


0.33 per year 


assuming 3 years life-span for untreated diseased hosts [81] 


b 


rate of transmission 


10 per year 


[82-84] 


w 


rate of regression 


0.25 per year 


return from active TB stage to latency [81] 


f 


fraction of TB Infections that lead to disease via fast 
progression 


0.1 


[11,46,50,84] 


k 


reduction of fast progression rate upon reinfection 
of latent hosts 


0.5 


[45,46,49] 


a 


rate of TB activation in latent hosts 


varied 





Initial conditions of model variables and values of model parameters. These values are chosen for all simulations unless indicated otherwise. 
doi:1 0.1 371/journal.pone.01 05721 .tool 
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Figure 2. Persistence as function of activation. A) Persistence, P, as a function of activation rate. B) Latent and infectious hosts, Land/, and the 
fraction of activation, a, as a function of activation rate. The left axis applies to L and /, the right axis to a. C) Persistence as function of fraction of hosts 
that activate. Also indicated in the figures is the optimal level of persistence, Pq, and the values for activation rate and fraction of activators at which 
this optimum occurs, ao and ao- All parameters are as given in table 1. 
doi:1 0.1 371/journal.pone.01 05721 .g002 



optimum occurs. Figure 3 shows the distribution of those values. 
Of most interest, for almost all parameter samples the fraction of 
hosts that activate remains well above the 10% found for MTB in 
the field. 

In text S 1 , we show further results by exploring how changes in 
each model parameter individually affect optimal persistence and 
activation. As expected, increased population size (though 
increased carrying capacity parameter or increased birth rate) 
leads to improved persistence, while decreased population size 
(through higher death rates) reduces persistence. Other model 
parameters have less impact on changes in persistence. The 
optimal fraction of activators as any of these parameters are 
changed individually is between 40%-70%, again above the value 
observed experimentally. 

We also explore in text SI how changes in the model structure, 
specifically the assumption of exponentially distributed natural 
lifespans, affect our results. We find that it shifts the persistence 
curve shown in fignre 2 slighdy, without affecting the overall 
results and still giving an optimal fraction of activators around 
50%. 

The deterministic persistence measure is well 
reproduced with a stochastic model 

Our persistence measure, P, is derived from a deterministic 
model. Of course, non-persistence, i.e. extinction, is an inherentiy 
stochastic process. While it is generally well known that larger 
population sizes lead to less chance of extinction [59-61], it is 
useful to directly test our deterministic measure with a stochastic 



model. To that end, we implemented the differential equations as 
a compartmental stochastic model, with transition rates of the 
deterministic model becoming transition probabilities [25,68]. We 
simulated the stochastic model using an efficient form of the 
Gillespie algorithm (the adaptive-T leap method as implemented in 
the R package adaptivetau [69,70]). Starting at the deterministic 
equilibrium state, we simulated the stochastic model for a fixed 
number of years and record the fraction of simulations for which 
at least one infectious or latent individual was still present at the 
end of the simulation. For the stochastic model, there is no 
discounting of the latent hosts by a factor a. Instead, persistence or 
extinction is a binary event, with extinction defined as no more 
latent or diseased infected hosts present and persistence if at least 
one of these hosts was still present. Figure 4 shows that despite this 
difference between P and the stochastic simulation, the fraction of 
simulations for which persistence was found in the stochastic 
model has a very similar functional shape as our deterministic 
persistence measure, P, with the optimum, Pq, occurring at more 
or less the same rate of activation. 

In Text SI, we use the stochastic model to investigate how 
population growth or decline influence P. We find that persistence 
improves in the presence of a growing population and worsens if 
the population dechnes. The shape of the persistence curve does 
not change in any important way. 

Persistence during epidemic cycles 

We have so far focused on MTB persistence at the endemic 
state. Equally important for pathogens is the ability to persist after 




Figure 3. Uncertainty analysis of model results. Impact of model parameter variations on A) optimal persistence, Po, B) optimal rate of 
activation, Qq, and C) optimal fraction of activation, zq. Model parameters were sampled in a range of 0.5-2 times the base parameters listed in 
table 1 . Sampling was done using a latin hypercube approach with 1 000 samples. Samples for which parameter combinations lead to a biologically 
unreasonable scenario, specifically natural mortality rate above birth rate, were discarded. For each sample, persistence as function of activation rate 
and activation fraction was determined (as shown in figure 2) and from this the optimal values for persistence and activation were obtained. 
doi:1 0.1 371/journal.pone.01 05721 .g003 
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introduction to a newly susceptible population. During TB's 
evolutionary history, it likely got introduced and re-introduced 
repeatedly into small susceptible subpopulations (e.g. new tribes/ 
villages). In the evolutionary context, persistence upon introduc- 
tion into a fiiUy susceptible population might therefore have played 
an important role. Nowadays, many areas of the world where TB 
is very rare consist of largely susceptible individuals - though for 
those groups the persistence idea explored here is likely not too 
important. Upon entering a fully susceptible population, patho- 
gens usually cause an epidemic outbreak, depleting the number of 
susceptibles. The outbreak is often followed by a fade-out of the 
disease once most susceptibles have been depleted. Extinction 



often occurs during this fade-out. For the pathogen to not go 
extinct, it needs to persist long enough until the number of 
susceptibles has built up again, usually leading to consecutive 
smaller outbreaks until the endemic state is reached [71,72]. We 
can quantify persistence during epidemic cycles by evaluating our 
expression for persistence not at the steady state but instead at the 
overall minimum occurring between introduction of the disease in 
a susceptible population and eventual attainment of the endemic 
equilibrium, i.e. we determine the overall minimum 
P,„= mm,[P(t)]. Figures shows the distribution of optimal 
persistence, optimal activation rate and optimal fraction of 
activating hosts for P„, using the same parameter sampling 
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Figure 4. Comparison of persistence with results from a stochastic model. The solid line shows persistence, P, as determined from the 
deterministic model, symbols show result from a stochastic model. For the stochastic simulation, we started at the deterministic steady state, 
simulated the model 10,000 times for 1000 years each and counted the numbers of simulations for which at least one latent or infectious host was 
still present after 1000 years. A population of size 100 was used, all parameter values are as reported in table 1. Note that because the absolute 
magnitude of P is arbitrary and scales with population size, N^, to allow comparison with the stochastic model we rescaled P to be between 0 and 1 . 
doi:1 0.1 371/journal.pone.01 05721. g004 
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Figure 5. Optimal persistence, P„„ activation rate and fraction for the non-steady state model. Everything as explained in caption for 
figure 3. 

doi:1 0.1 371 /journal.pone.01 05721 .g005 



approach as for figure 3. Comparing the results with those found 
earher for P at the steady state (figure 3), one sees that the results 
are very similar. The main reason for the similarity is that TB has 
a relatively "slow" disease dynamics [73], without pronounced 
outbreak peaks and minima. Therefore, for most parameter 
values, the disease reaches the steady state without a large 
contraction after the first outbreak, leading to essentially the same 
results as for the steady state. 

Again, as for P above, we show additional results for P„ as 
function of individual parameters in Text SI. The results are 
almost identical to those obtained for P, for the reason just 
explained. 

Discussion 

Recent in-depth studies of MTB genetic sequences have shown 
a wide diversity between strains [17]. While harder to determine, 
there is also accumulating evidence that this genetic diversity 
results in phenotypic diversity [17,74], suggesting that MTB 
evolution is shaped by local selection pressures. A recent study 
indicates that disease activation rate differs between lineages, 
suggesting that this phenotype is under evolutionary selection [75]. 
We used a mathematical model to investigate the role of activation 
rate and latency duration on the ability of MTB to persist in a host 
population. Our results support the previously proposed idea that 
the prolonged latency observed for TB infections might provide 
MTB an evolutionary advantage, namely improved persistence in 
a host population [23,41]. We found that an intermediate rate of 
disease activation is optimal for persistence. 

Interestingly, our model suggests that for optimal persistence, 
the fraction of hosts that eventually become cases is in the 20%- 
80% range with approximately 50% as the most likely value. This 
range is higher than the —10% generally seen for TB, suggesting 
that the host immune response plays some role in keeping TB 
disease in check, lowering activation rate below what would be the 
evolutionary optimal level for the pathogen. It is likely that some 
level of co-evolution between pathogen and host occurred and that 
humans who have been exposed to MTB for a long time evolved 
some level of resistance that potentially prevents MTB from 
reaching its evolutionary optimal activation rate [12]. This also 
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